#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "funciones.h"

struct matrix {
		int n;
		double **a;
		double *x;
		double *y;
		double *b;
};

struct matrix *create_matrix(unsigned int n){
	int i = 0;
	struct matrix *m = NULL;
	assert(n > 0);
	
	m = (struct matrix *)calloc(1, sizeof(struct matrix));
	m->n = n;
	
	m->a = (double **)calloc(n, sizeof(double));
	for(i = 0 ; i < n ; i++){
			m->a[i] = (double *)calloc(n, sizeof(double));
	}
	m->x = (double *)calloc(n, sizeof(double));
	m->y = (double *)calloc(n, sizeof(double));
	m->b = (double *)calloc(n, sizeof(double));
	
	return m;
}

struct matrix *destruct_matrix(struct matrix *m){
	int i = 0;
	
	assert(m != NULL);
	
	for(i = 0 ; i < m->n ; i++){
		free(m->a[i]);
		m->a[i] = NULL;
	}
	free(m->a);
	m->a = NULL;
	free(m->x);
	m->x = NULL;
	free(m->y);
	m->y = NULL;
	free(m->b);
	m->b = NULL;
	free(m);
	m = NULL;

	return m;	
}

double distancia (struct matrix *m) {
	int i = 0;
	double d = 0.0;
	while (i < m->n) {
		d = d + sqrt((m->x[i] - m->y[i])*(m->x[i] - m->y[i]));
		i = i+1;
		}
return d;
}

void iteracion (struct matrix * m) {
	int i, j;
	double sigma = 0.0;
	
	for (i = 0; i < m->n; i++) {
		sigma = 0.0;
		for (j = 0; j < m->n; j++) {
			if (j != i) {
				sigma = (sigma + (m->a[i][j] * m->x[j]));
			}
		}
	m->y[i] = (m->b[i]-sigma)/m->a[i][i];
	}
}

void jacobi (struct matrix *m, double delta, int M, double *norma) {
	int i, w = 1;
	iteracion (m);
	*norma = distancia(m);
	while (w <= M && *norma >= delta) {
		for (i = 0; i < m->n; i++) {
		m->x[i] = m->y[i];
		}
		iteracion (m);
		printf("El vector de la %d° aproximacion es: \n", w);
		for (i = 0; i < m->n; i++) {
		printf("x[%d] = %.10lf\n",i,m->x[i]);
		}
		w = w + 1;
		*norma = distancia (m);
		printf("La distancia entre las dos ultimas aproximaciones es: %.10lf\n\n",*norma);
	}
}

struct matrix *matrices (void) {
	
int i, j, n;
struct matrix *m = NULL;
char linea[1024];
char *r;
FILE *datos;

        datos = fopen("datos","r");
        if ( datos == NULL ) {
                printf ("No se puede abrir el archivo\n");
                exit(1);
        }
        
        r = fgets(linea, 1024, datos);
        n = atoi(r);

		m = create_matrix(n);
        for (i=0; i<n; i++) {
          for (j=0; j<n; j++) {
            r = fgets(linea, 1024, datos);
            if (i == j && !strncmp(r,"0",1)) {
				printf("El metodo de Jacobi no es aplicable a la matriz. \n");
				exit(0);
			}
            m->a[i][j]= atof(r);
          }
        }
             
        for (i=0; i<n; i++) {
          r = fgets(linea, 1024, datos);
          m->x[i]= atof(r);
        }
        for (i=0; i<n; i++) {
          r = fgets(linea, 1024, datos);
          m->b[i]= atof(r);
        }
        
        fclose(datos);
        
        printf("Dimension: %d\n", n);
        printf("La matriz A ingresada es: \n");
        for (i=0; i<n; i++) {
          for (j=0; j<n; j++) {
            printf ("%.2f ",m->a[i][j]);
          }
          printf("\n");
        }
        printf("El vector solucion inicial es: \n");
        for (i=0; i<n; i++) {
          printf ("%f\n",m->x[i]);
        }
        printf("El vector b es: \n");
        for (i=0; i<n; i++) {
          printf ("%f\n",m->b[i]);
        }
		
return m;
}

int jacobifinal (void) {
	int M;
	double norma = 0.0;
	double delta;
	struct matrix *m =  NULL;
		
    m = matrices();
	
	printf("Ingrese el error permitido delta: \n");
    scanf("%lf",&delta);
        
    printf("Ingrese el maximo de iteraciones permitidas: \n");
    scanf("%d",&M);
        
	iteracion(m);
	jacobi(m, delta,M,&norma);
	
	m = destruct_matrix(m);
	
return 0;
}

void iteracion_gs (struct matrix * m) {
	int i, j;
	double sigma1 = 0.0;
	double sigma2 = 0.0;
	
	for (i = 0; i < m->n; i++) {
		sigma1 = 0.0;
		sigma2 = 0.0;
		for (j = 0; j < i; j++) {
			sigma1 = (sigma1 + (m->a[i][j] * m->y[j]));
			}
		for (j = i+1; j < m->n; j++) {
			sigma2 = (sigma2 + (m->a[i][j] * m->x[j]));
			}
	m->y[i] = (m->b[i]-sigma1-sigma2)/m->a[i][i];
	}
}

void gauss_seidel (struct matrix *m, double delta, int M, double *norma) {
	int i, w = 1;
	iteracion_gs (m);
	*norma = distancia(m);
	while (w <= M && *norma >= delta) {
		for (i = 0; i < m->n; i++) {
		m->x[i] = m->y[i];
		}
		iteracion_gs (m);
		printf("El vector de la %d° aproximacion es: \n", w);
		for (i = 0; i < m->n; i++) {
		printf("x[%d] = %.10lf\n",i,m->x[i]);
		}
		w = w + 1;
		*norma = distancia (m);
		printf("La distancia entre las dos ultimas aproximaciones es: %.10lf\n\n",*norma);
	}
}

int gauss_seidel_final (void) {
	
	int M;
	double norma = 0.0;
	double delta;
	struct matrix *m =  NULL;
		
    m = matrices();
	
	printf("Ingrese el error permitido delta: \n");
    scanf("%lf",&delta);
        
    printf("Ingrese el maximo de iteraciones permitidas: \n");
    scanf("%d",&M);
        
	iteracion_gs(m); 
	gauss_seidel(m,delta,M,&norma);
	
	m = destruct_matrix(m);
	
return 0;
}
